Host and Environmental drivers of microbiome: Beta diversity

This script carries out the beta diversity analyses (community composition and structure) for the 2016 Cork great tit microbiome data.

Overall plan" * model 1: all variables of interest, nest as blocking factor * model 2: all variables of interest but age-habitat interaction, nest blocking * model 3: adult samples, all variables of interest, nest blocking * then check are results different when individuals w/ 2 samples are dropped

Gabrielle says not to remove individuals sampled twice

Set-up and libraries

library(devtools)

Loading required package: usethis

Attaching package: ‘devtools’

The following object is masked from ‘package:permute’:

    check

install_github("pmartinezarbizu/pairwiseAdonis/pairwiseAdonis")

Skipping install of 'pairwiseAdonis' from a github remote, the SHA1 (ece560d2) has not changed since last install.
  Use `force = TRUE` to force installation

phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 54343 taxa and 246 samples ]
sample_data() Sample Data:       [ 246 samples by 50 sample variables ]
tax_table()   Taxonomy Table:    [ 54343 taxa by 7 taxonomic ranks ]
refseq()      DNAStringSet:      [ 54343 reference sequences ]

Data cleaning

Remove individuals measured twice * set.seed so same individuals dropped * no adults measured twice

#n_occur <- data.frame(table(metadata.NoDups$bird.ID))
#n_occur[n_occur$Freq > 1,]

set.seed(1189)
metadata.NoDups <- metadata %>% 
    group_by(bird.ID) %>%
    sample_n(1)

dropped.samples <- setdiff(metadata$BIOM.ID, metadata.NoDups$BIOM.ID)

# metadata.even <- metadata %>% 
#     group_by(bird.ID) %>%
#     sample_n(min(table(metadata$bird.ID)))

table(metadata$ageBinned)

1week 2week adult 
   81   114    51 

table(metadata.NoDups$ageBinned)

1week 2week adult 
   55    91    51 

Do i need to scale numeric variables?

numeric.predictors <- c("QubitDNA","Tarsus","Weight", "wing", "broodSizeWhenSampled", "broodSizeMax", "totalFledge", "clutchSize", "layDateFirst", "numberDeadPreRinged", "numberDeadPostRinged","scaled.mass","scaled.mass.wing.adult", "scaled.mass.tarsus.adult","scaled.mass.chick", "scaled.mass.tarsus", "scaled.mass.wing", "DistanceToEdge")

# centre and scale numeric variables ie. subtract mean and divide by st. deviation
metadata.NoDups.scaled <- metadata.NoDups

Error: object 'metadata.NoDups' not found

All birds: CoDa, w/ duplicates

Compositional method, use clr transform

Keeping repeated measures here but using blocking factor to control for repeated samples.

# Filter rare taxa
phylo.knowles <- filter_taxa(phylo.spring, function(x) sum(x > 1) > (0.05*length(x)), TRUE)

# phylo.knowles.clr <- phylo.knowles
# phylo.knowles.clr@otu_table <- otu_table(clr(phylo.knowles@otu_table), taxa_are_rows = F)
# aitchison.dist.dups <- phyloseq::distance(phylo.knowles.clr, method = "euclidean")
pk.otu.clr.dups <- clr(phylo.knowles@otu_table)
aitchison.dist.dups <- vegdist(pk.otu.clr.dups, method = "euclid")

# phylo.TSS <- transform_sample_counts(phylo.knowles, function(x) x/sum(x)) # normalise read counts w/ Total-Sum Scaling
# 
# BCdist <- phyloseq::distance(phylo.TSS, method="bray")

All variables except for habitat have heterogenous dispersions

## H0= No difference in dispersion between groups
# calc dispersion, using distance measure

variables <- c("ageBinned", "habitat", "layDateFirst", "broodSizeWhenSampled", "DistanceToEdge", "SequencePlate")

for(i in variables){ # works
  dispersion <- betadisper(aitchison.dist.dups, metadata.scaled[,i])  #, bias.adjust = T
  print(i) # print variable being tested
  print(permutest(dispersion, pairwise=FALSE, permutations=1000))
  cat("\n") # print line break, makes it easier to read
}

#hist(metadata$layDateFirst)
#ist(sqrt(metadata$DistanceToEdge))

dispersion1 <- betadisper(aitchison.dist.dups, metadata$ageBinned)  #, bias.adjust = T
print(permutest(dispersion1, pairwise=FALSE, permutations=1000))

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 1000

Response: Distances
           Df Sum Sq Mean Sq      F N.Perm  Pr(>F)  
Groups      2  107.0    53.5 3.7152   1000 0.02498 *
Residuals 243 3499.2    14.4                        
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

dispersion2 <- betadisper(aitchison.dist.dups, metadata$habitat)  #, bias.adjust = T
print(permutest(dispersion2, pairwise=FALSE, permutations=1000))

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 1000

Response: Distances
           Df Sum Sq Mean Sq      F N.Perm Pr(>F)
Groups      1    5.7  5.7331 0.3846   1000 0.5504
Residuals 244 3637.6 14.9082                     

dispersion3 <- betadisper(aitchison.dist.dups, (metadata$layDateFirst))  #, bias.adjust = T
print(permutest(dispersion3, pairwise=FALSE, permutations=1000))

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 1000

Response: Distances
           Df Sum Sq Mean Sq      F N.Perm   Pr(>F)    
Groups     27 1047.2  38.785 2.6146   1000 0.000999 ***
Residuals 218 3233.8  14.834                           
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

dispersion4 <- betadisper(aitchison.dist.dups, metadata$broodSizeWhenSampled)  #, bias.adjust = T
print(permutest(dispersion4, pairwise=FALSE, permutations=1000))

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 1000

Response: Distances
           Df Sum Sq Mean Sq     F N.Perm  Pr(>F)  
Groups      6  220.0  36.674 2.505   1000 0.02498 *
Residuals 239 3499.1  14.641                       
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

dispersion5 <- betadisper(aitchison.dist.dups, (metadata$DistanceToEdge))  #, bias.adjust = T
print(permutest(dispersion5, pairwise=FALSE, permutations=1000))

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 1000

Response: Distances
           Df Sum Sq Mean Sq      F N.Perm   Pr(>F)    
Groups     60 2515.5  41.924 4.0906   1000 0.000999 ***
Residuals 185 1896.0  10.249                           
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

dispersion6 <- betadisper(aitchison.dist.dups, metadata$SequencePlate)  #, bias.adjust = T
print(permutest(dispersion6, pairwise=FALSE, permutations=1000))

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 1000

Response: Distances
           Df  Sum Sq Mean Sq      F N.Perm   Pr(>F)    
Groups      4  541.68 135.419 10.509   1000 0.000999 ***
Residuals 241 3105.49  12.886                           
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

TukeyHSD(dispersion1)

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = distances ~ group, data = df)

$group
                  diff       lwr        upr     p adj
2week-1week -1.4854366 -2.785847 -0.1850262 0.0205827
adult-1week -1.1177060 -2.717328  0.4819163 0.2277401
adult-2week  0.3677306 -1.139788  1.8752487 0.8334642

TukeyHSD(dispersion2)

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = distances ~ group, data = df)

$group
                       diff        lwr      upr     p adj
deciduous-conifer 0.3765808 -0.8195671 1.572729 0.5357527

#TukeyHSD(dispersion3)
#TukeyHSD(dispersion4)
#TukeyHSD(dispersion5)

Plot dispersions

perms.dups <- with(metadata, how(nperm = 1000, blocks = nest))

all.adonis.dups.fixed <- adonis2(aitchison.dist.dups ~ ageBinned + habitat + layDateFirst + broodSizeWhenSampled + DistanceToEdge + SequencePlate, by="margin", method="euclidian", data = metadata.scaled, permutations = perms.dups)

pairwiseAdonis::pairwise.adonis2(aitchison.dist.dups ~ ageBinned + habitat, data = metadata.scaled)

$parent_call
[1] "aitchison.dist.dups ~ ageBinned + habitat , strata = Null , permutations 999"

$`2week_vs_adult`
           Df SumOfSqs      R2      F Pr(>F)    
ageBinned   1      647 0.01699 2.8293  0.001 ***
habitat     1      395 0.01036 1.7263  0.001 ***
Residual  162    37034 0.97265                  
Total     164    38075 1.00000                  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

$`2week_vs_1week`
           Df SumOfSqs      R2      F Pr(>F)    
ageBinned   1      635 0.01317 2.5796  0.001 ***
habitat     1      335 0.00695 1.3625  0.027 *  
Residual  192    47228 0.97988                  
Total     194    48198 1.00000                  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

$adult_vs_1week
           Df SumOfSqs      R2      F Pr(>F)    
ageBinned   1      563 0.01645 2.1756  0.001 ***
habitat     1      270 0.00788 1.0425  0.365    
Residual  129    33369 0.97566                  
Total     131    34201 1.00000                  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

attr(,"class")
[1] "pwadstrata" "list"      

#?pairwise.adonis2()

Try same as above but strata by sequence plate and then inlcude individual ID

# perms.dups <- with(metadata, how(nperm = 1000, blocks = SequencePlate))
# 
# all.adonis.dups.fixed <- adonis2(aitchison.dist.dups ~ ageBinned + habitat + layDateFirst + broodSizeWhenSampled + DistanceToEdge + bird.ID, by="margin", method="euclidian", data = metadata.scaled, permutations = perms.dups)
# 
# all.adonis.dups.fixed

Including bird.ID as a fixed effect and blocking by nest suggests age is in fact a significant factor though accounts for only 0.7% of variation, while bird id accounts for 75%, but is non-significant. Sequence plate is significant, accounting for 2.5%, brood size is not significant while habitat, lay date and distance to edge could not be estimated.

# perms.dups <- with(metadata, how(nperm = 1000, blocks = nest))
# 
# all.adonis.dups.fixed <- adonis2(aitchison.dist.dups ~ ageBinned + habitat + layDateFirst + broodSizeWhenSampled + DistanceToEdge + SequencePlate + bird.ID, by="margin", method="euclidian", data = metadata.scaled, permutations = perms.dups)
# 
# all.adonis.dups.fixed

# perms.dups <- with(metadata, how(nperm = 1000, blocks = bird.ID))
# 
# all.adonis.dups.fixed <- adonis2(aitchison.dist.dups ~ ageBinned + habitat + layDateFirst + broodSizeWhenSampled + DistanceToEdge + SequencePlate, by="margin", method="euclidian", data = metadata.scaled, permutations = perms.dups)
# 
# all.adonis.dups.fixed

perms.dups <- with(metadata, how(nperm = 1000, blocks = nest))

all.adonis.dups.int <- adonis2(aitchison.dist.dups ~ ageBinned*habitat + layDateFirst + broodSizeWhenSampled + DistanceToEdge + SequencePlate, by="margin", data = metadata.scaled, permutations = perms.dups) #
all.adonis.dups.int

Permutation test for adonis under reduced model
Marginal effects of terms
Blocks:  nest 
Permutation: free
Number of permutations: 1000

adonis2(formula = aitchison.dist.dups ~ ageBinned * habitat + layDateFirst + broodSizeWhenSampled + DistanceToEdge + SequencePlate, data = metadata.scaled, permutations = perms.dups, by = "margin")
                      Df SumOfSqs      R2      F   Pr(>F)    
layDateFirst           1      420 0.00693 1.8300 0.089910 .  
broodSizeWhenSampled   1      305 0.00504 1.3312 0.016983 *  
DistanceToEdge         1      502 0.00830 2.1901 0.000999 ***
SequencePlate          4     3634 0.06002 3.9609 0.000999 ***
ageBinned:habitat      2      595 0.00982 1.2960 0.010989 *  
Residual             233    53444 0.88259                    
Total                245    60554 1.00000                    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Interaction calc by term rather than margin to get estimates for main effects as well as interaction term

perms.dups <- with(metadata, how(nperm = 1000, blocks = nest))

adonis2(aitchison.dist.dups ~ ageBinned*habitat + layDateFirst + broodSizeWhenSampled + DistanceToEdge + SequencePlate, by="term", data = metadata.scaled, permutations = perms.dups) #

Permutation test for adonis under reduced model
Terms added sequentially (first to last)
Blocks:  nest 
Permutation: free
Number of permutations: 1000

adonis2(formula = aitchison.dist.dups ~ ageBinned * habitat + layDateFirst + broodSizeWhenSampled + DistanceToEdge + SequencePlate, data = metadata.scaled, permutations = perms.dups, by = "term")
                      Df SumOfSqs      R2      F   Pr(>F)    
ageBinned              2     1239 0.02046 2.7003 0.000999 ***
habitat                1      347 0.00573 1.5124 0.000999 ***
layDateFirst           1      490 0.00810 2.1377 0.000999 ***
broodSizeWhenSampled   1      311 0.00514 1.3578 0.078921 .  
DistanceToEdge         1      475 0.00784 2.0689 0.000999 ***
SequencePlate          4     3653 0.06033 3.9816 0.000999 ***
ageBinned:habitat      2      595 0.00982 1.2960 0.009990 ** 
Residual             233    53444 0.88259                    
Total                245    60554 1.00000                    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

check influence of unbalanced design

table(metadata$ageBinned, metadata$habitat)

       
        conifer deciduous
  1week      16        65
  2week      27        87
  adult       8        43

For loop sample equal numbers of from each habitat and age combination

Plot PCA

aitch.pca <- prcomp(aitchison.dist.dups)
#plot(aitch.pca)
biplot(aitch.pca)

autoplot(aitch.pca, data = metadata)#, colour = "ageBinned", loadings = T, loadings.label = T, scale = 0)

Code from: https://huboqiang.cn/2016/03/03/RscatterPlotPCA

df_out <- as.data.frame(aitch.pca$x)
df_out$ageBinned <- metadata$ageBinned
df_out$habitat <- metadata$habitat
df_out$layDateFirst <- metadata$layDateFirst
df_out$broodSizeWhenSampled <- metadata$broodSizeWhenSampled

head(df_out)

Calculate the PC loadings for labels * took these from different sites and they disagree alot

# Sum the total variance
d.mvar <- sum(aitch.pca$sdev^2)
# Calculate the PC1 and PC2 variance
PC1.label <- paste("PC1: ", round(sum(aitch.pca$sdev[1]^2)/d.mvar, 3)*100,"%")
PC2.label <- paste("PC2: ", round(sum(aitch.pca$sdev[2]^2)/d.mvar, 3)*100,"%")

#percentage <- round(aitch.pca$sdev / sum(aitch.pca$sdev) * 100, 2)
#percentage <- paste( colnames(df_out), "(", paste( as.character(percentage), "%", ")", sep="") )

#screeplot(aitch.pca)

Plot * Use same name and label for fill and shape in order to overlap legends * ellipse default is 95% confidence level for t distribution * can draw euclid ellipse where level = radius, but how to choose radius? * alpha specifies the transparency * https://stats.stackexchange.com/questions/217374/real-meaning-of-confidence-ellipse

p <- ggplot(df_out,aes(x=PC1,y=PC2, color=ageBinned ))
p <- p + geom_point(aes(shape = habitat)) + xlab(PC1.label) + ylab(PC2.label) #, show.legend = F
p <- p + stat_ellipse(geom = "polygon", type="t", alpha=0.1, 
                      aes(group = interaction(ageBinned, habitat), fill = habitat, colour = ageBinned))
p <- p + scale_fill_discrete(name = "Habitat", labels = c("Conifer", "Deciduous"))
p <- p + scale_shape_discrete(name = "Habitat", labels = c("Conifer", "Deciduous"))
p <- p + scale_color_discrete(name = "Age", labels = c("Day 8", "Day 15", "Adult"))
#p <- p + scale_shape_discrete(element_blank())#, labels = NULL)
#p + facet_grid(habitat~.)

#theme<-theme(panel.background = element_blank(), panel.border=element_rect(fill=NA), panel.grid.major = element_blank(), panel.grid.minor = element_blank(),strip.background=element_blank(), axis.text.x=element_text(colour="black"), axis.text.y=element_text(colour="black"), axis.ticks=element_line(colour="black"), plot.margin=unit(c(1,1,1,1),"line"))

p+theme_classic()

Try euclidian ellipses, specify radius with level

p <- ggplot(df_out,aes(x=PC1,y=PC2, color=ageBinned ))
p <- p + geom_point(aes(shape = habitat)) + xlab(PC1.label) + ylab(PC2.label) #, show.legend = F
p <- p + stat_ellipse(geom = "polygon", type="euclid", alpha=0.1,level = 10,
                      aes(group = interaction(ageBinned, habitat), fill = habitat, colour = ageBinned))
p <- p + scale_fill_discrete(name = "Habitat", labels = c("Conifer", "Deciduous"))
p <- p + scale_shape_discrete(name = "Habitat", labels = c("Conifer", "Deciduous"))
p <- p + scale_color_discrete(name = "Age", labels = c("Day 8", "Day 15", "Adult"))
#p <- p + scale_shape_discrete(element_blank())#, labels = NULL)
#p + facet_grid(habitat~.)

#theme<-theme(panel.background = element_blank(), panel.border=element_rect(fill=NA), panel.grid.major = element_blank(), panel.grid.minor = element_blank(),strip.background=element_blank(), axis.text.x=element_text(colour="black"), axis.text.y=element_text(colour="black"), axis.ticks=element_line(colour="black"), plot.margin=unit(c(1,1,1,1),"line"))

p+theme_classic()

facet.labels <- c("Conifer", "Deciduous")
names(facet.labels) <- c("conifer", "deciduous")
p + theme_classic() + facet_grid(habitat~., labeller = labeller(habitat = facet.labels))

# plot_ordination(phylo.TSS, ordBC, color = "ageBinned") + 
#   geom_point(size=3) + 
#   ggtitle("Ord: Bray-Curtis") + stat_ellipse(geom = "polygon", type="t", alpha=0.1, aes(fill=ageBinned))+ 
#   theme_bw() +  theme(panel.grid.major = element_blank(),
#                       panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"), 
#                       panel.border = element_rect(linetype = "solid", colour = "black", size=.8)) +
#                       theme(text=element_text(size=14,  family="serif"), axis.ticks = element_line(colour = "black", size = .7))+
#                       geom_point()

Plot features that contribute to classification

df_out_r <- as.data.frame(aitch.pca$rotation)
df_out_r$feature <- row.names(df_out_r)
df_out_r$ageBinned <- metadata$ageBinned
df_out_r$habitat <- metadata$habitat
df_out_r$layDateFirst <- metadata$layDateFirst
df_out_r$broodSizeWhenSampled <- metadata$broodSizeWhenSampled
df_out_r$DistanceToEdge <- metadata$DistanceToEdge


df_out_r

p<-ggplot(df_out_r,aes(x=PC1,y=PC2,label=ageBinned,color=habitat ))
#p<-p+geom_point()+theme + geom_text(size=3) + theme(legend.position = "none")
p

Plot one plot for each habitat ((ellipses each age) * change aplha to 0.05?

p <- ggplot(df_out,aes(x=PC1,y=PC2, color=ageBinned ))
p <- p + geom_point() + xlab(PC1.label) + ylab(PC2.label) #, show.legend = F
p <- p + stat_ellipse(geom = "polygon", type="t", alpha=0.1, 
                      aes(group = interaction(ageBinned), fill = ageBinned, colour = ageBinned))
p <- p + scale_color_discrete(name = "Age", labels = c("Day 8", "Day 15", "Adult"))
p <- p + scale_fill_discrete(element_blank(), labels = NULL)
#p + facet_grid(habitat~.)

p+theme_classic()

facet.labels <- c("Conifer", "Deciduous")
names(facet.labels) <- c("conifer", "deciduous")
p + theme_classic() + facet_grid(habitat~., labeller = labeller(habitat = facet.labels))

one plot for each age (ellipses each habitat)

p <- ggplot(df_out,aes(x=PC1,y=PC2, color=habitat ))
p <- p + geom_point() + xlab(PC1.label) + ylab(PC2.label) #, show.legend = F
p <- p + stat_ellipse(geom = "polygon", type="t", alpha=0.1, 
                      aes(group = (habitat), fill = habitat, colour = habitat))
p <- p + scale_color_discrete(name = "Habitat", labels = c("Conifer", "Deciduous"))
p <- p + scale_fill_discrete(element_blank(), labels = NULL)
#p + facet_grid(habitat~.)

p+theme_classic()

facet.labels <- c("Day 8", "Day 15", "Adult")
names(facet.labels) <- c("1week", "2week", "adult")
p + theme_classic() + facet_grid(ageBinned~., labeller = labeller(ageBinned = facet.labels))

JQ requested plots

ellipses for just age

one plot for each age (ellipses each habitat)

p <- ggplot(df_out,aes(x=PC1,y=PC2, color=ageBinned ))
p <- p + geom_point() + xlab(PC1.label) + ylab(PC2.label) #, show.legend = F
p <- p + stat_ellipse(geom = "polygon", type="t", alpha=0.1, 
                      aes(group = (ageBinned), fill = ageBinned, colour = ageBinned))
p <- p + scale_color_discrete(name = "Age", labels = c("D8", "D15", "Adult"))
p <- p + scale_fill_discrete(element_blank(), labels = NULL)
#p + facet_grid(habitat~.)

p+theme_classic()

ellipses for just habitat

p <- ggplot(df_out,aes(x=PC1,y=PC2, color=habitat ))
p <- p + geom_point() + xlab(PC1.label) + ylab(PC2.label) #, show.legend = F
p <- p + stat_ellipse(geom = "polygon", type="t", alpha=0.1, 
                      aes(group = (habitat), fill = habitat, colour = habitat))
p <- p + scale_color_discrete(name = "Habitat", labels = c("Conifer", "Deciduous"))
p <- p + scale_fill_discrete(element_blank(), labels = NULL)
#p + facet_grid(habitat~.)

p+theme_classic()

p.habitat <- p

Combine age and habitat plots

  theme_bw(base_size = base_size, base_family = base_family) %+replace%
    theme(
      legend.position = "none",
      panel.border     = element_blank(),
      axis.line        = element_line(colour = "black"),
      panel.grid.major=element_line(colour="grey", size=0.5, 3),
      panel.grid.major.x = element_blank(),
      panel.grid.major.y = element_blank(),
      panel.grid.minor = element_blank(),
      panel.grid.minor.x = element_blank(),
      panel.grid.minor.y = element_blank(),
      strip.background = element_blank(),
      legend.key       = element_blank(),
      plot.title = element_text(hjust = 0.5, face = "bold") # centre and bold title
    )

Error in theme_grey(base_size = base_size, base_family = base_family,  : 
  object 'base_size' not found

brood size

envfit1 <- envfit(aitch.pca,metadata$broodSizeWhenSampled, main="", col="red")
#plot(aitch.pca, type = "n")
plot(envfit1)

Error in strwidth(labels, ...) : plot.new has not been called yet

Smooth surface/surf plot

# from Jennys email + https://chrischizinski.github.io/rstats/ordisurf/
species.scores <- as.data.frame(scores(aitch.pca, "species")) # have to keep species here as the taxa

Error in names(x) <- paste("V", seq_along(x), sep = "") : 
  'names' attribute [1] must be the same length as the vector [0]

https://userweb.eng.gla.ac.uk/umer.ijaz/bioinformatics/ecological.html

#p<-p+theme(legend.key = element_blank(),  #removes the box around each legend item
#              legend.position = "bottom", #legend at the bottom
#              legend.direction = "horizontal",
#              legend.box = "horizontal",
#              legend.box.just = "centre")
p

Error in rep(x[[i]], n) : 
  attempt to replicate an object of type 'closure'

Adult birds: CoDa

Aitchison distance is supposedly robust to subsetting- can i therefore just subset the adult birds from previous distance calculation?

phylo.kn.adult <- filter_taxa(phylo.adults, function(x) sum(x > 1) > (0.05*length(x)), TRUE)

# phylo.TSS.adults <- transform_sample_counts(phylo.kn.adult, function(x) x/sum(x))
# 
# BCdist.adults <- phyloseq::distance(phylo.TSS.adults, method="bray")
# JDdist.adults <- phyloseq::distance(phylo.TSS.adults, method="jaccard")

pk.otu.clr.adult <- clr(phylo.kn.adult@otu_table)
aitchison.dist.adult <- vegdist(pk.otu.clr.adult, method = "euclidean")

Check dispersion

Homogenous dispersion: sex, habitat, distance is almost non-homo Non-homogenous: ageDays, …all others

## H0= No difference in dispersion between groups
# calc dispersion, using distance measure

variables <- c("ageDays", "Sex", "habitat", "layDateFirst", "broodSizeWhenSampled", "DistanceToEdge", "SequencePlate")

for(i in variables){ # works
  dispersion <- betadisper(aitchison.dist.adult, adults.meta.scaled[,i])
  print(i) # print variable being tested
  print(permutest(dispersion, pairwise=FALSE, permutations=1000))
  cat("\n") # print line break, makes it easier to read
}

Modelling: adult aitchison

Sex * habitat * broodSizeWhenSampled + Sex * habitat * layDateFirst + habitat * DistanceToEdge + (1|nest) + (1|SequencePlate)

perms.adult <- with(adults.meta, how(nperm = 1000, blocks = nest))

adult.adonis <- adonis2(aitchison.dist.adult ~ ageDays + Sex + habitat + DistanceToEdge + layDateFirst + broodSizeWhenSampled + SequencePlate, by = "margin", data = adults.meta.scaled, permutations = perms.adult)
adult.adonis

Permutation test for adonis under reduced model
Marginal effects of terms
Blocks:  nest 
Permutation: free
Number of permutations: 1000

adonis2(formula = aitchison.dist.adult ~ ageDays + Sex + habitat + DistanceToEdge + layDateFirst + broodSizeWhenSampled + SequencePlate, data = adults.meta.scaled, permutations = perms.adult, by = "margin")
                     Df SumOfSqs      R2      F   Pr(>F)   
ageDays               1    271.0 0.01968 1.0475 0.629371   
Sex                   1    213.0 0.01546 0.8231 0.844156   
habitat               1    290.3 0.02108 1.1221 0.449550   
DistanceToEdge        1    327.0 0.02375 1.2640 0.315684   
layDateFirst          1    323.6 0.02350 1.2507 0.145854   
broodSizeWhenSampled  1    223.3 0.01621 0.8631 0.001998 **
SequencePlate         3   1462.1 0.10616 1.8836 0.003996 **
Residual             41  10608.1 0.77026                   
Total                50  13772.1 1.00000                   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#adonis2(JDdist.adults ~ Sex + habitat + DistanceToEdge + layDateFirst, by = "margin", data = adults.meta, permutations = perms.adult)

Results

Table of all results for fixed and interaction PERMANOVA plus adult

#make df of lmer output
all.bdiv <- as_tibble(rbind(all.adonis.dups.fixed, all.adonis.dups.int, adult.adonis), rownames="Independent variables") %>%
  dplyr::rename("P_estimate"="Pr(>F)") %>% 
  mutate_if(is.numeric, round, 3) %>%
  mutate(P_estimate=ifelse(P_estimate==0,"<0.001",P_estimate)) %>% 
  mutate(P_estimate=ifelse(P_estimate<=0.05,str_c(P_estimate," *"),P_estimate)) %>%
  mutate(P_estimate=ifelse(P_estimate>=0.05 & P_estimate<=0.06,str_c(P_estimate,"  ."),P_estimate))


# make df into kable
## this is in html, dosnt render in word doc # but can copy-paste
all.bdiv.results <- kable(all.bdiv, format = "html", table.attr = "style = \"color: black;\"") %>%
  kableExtra::kable_styling(full_width = F) #%>%
 # kableExtra::group_rows("(a) Top model, proteobacteria",1,nrow(proteobacteria.modelAvg.df)) %>%
  #
  #save_kable("alpha-kable__________.png") # this line saves as .png in Reports/

all.bdiv.results